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ABSTRACT 

The statistics of dark matter halos is an essential component of precision cosmology. The mass distribution 
of halos, as specified by the halo mass function, is a key input for several cosmological probes. The sizes of 
TV-body simulations are now such that, for the most part, results need no longer be statistics-limited, but are still 
subject to various systematic uncertainties. Discrepancies in the results of simulation campaigns for the halo 
mass function remain in excess of statistical uncertainties and of roughly the same size as the error limits set 
by near-future observations; we investigate and discuss some of the reasons for these differences. Quantifying 
error sources and compensating for them as appropriate, we carry out a high-statistics study of dark matter halos 
from 67 TV-body simulations to investigate the mass function and its evolution for a reference ACDM cosmology 
and for a set of wCDM cosmologies. For the reference ACDM cosmology (close to WMAP5), we quantify the 
breaking of universality in the form of the mass function as a function of redshift, finding an evolution of as 
much as 10% away from the universal form between redshifts z = and z = 2. For cosmologies very close to this 
reference we provide a fitting formula to our results for the (evolving) ACDM mass function over a mass range 
of 6 ■ 10 1 ' - 3 • 10 15 M Q to an estimated accuracy of about 2%. The set of wCDM cosmologies is taken from the 
Coyote Universe simulation suite. The mass functions from this suite (which includes a ACDM cosmology and 
others with w ~ -1) are described by the fitting formula for the reference ACDM case at an accuracy level of 10%, 
but with clear systematic deviations. We argue that, as a consequence, fitting formulae based on a universal form 
for the mass function may have limited utility in high precision cosmological applications. 

Subject headings: Cosmology: large-scale structure of universe — methods: TV-body simulations 



1. INTRODUCTION 

The current paradigm for the formation of cosmological struc- 
ture is based on the gravitational amplification of primordial 
density fluctuations in an expanding Universe. The nonlinear 
transformation of dark matter overdensities - via a hierarchi- 
cal dynamical process - into clumpy distributions called halos, 
and the subsequent infall of baryons leading to the formation 
of stars and galaxies within these halos, rounds out the present 
picture of the formation of observed structure. Although there 
is no precise mathematical definition of a 'halo', several opera- 
tional definitions - depending on the particular applications of 
interest - have been employed in practice. 

The spatio-temporal statistics of halos and sub-halos, as well 
as of their mass distribution (and its evolution), together pro- 
vide most of the descriptive framework within which fit all of 
the structure formation-based probes of cosmology. The mass 
function alone is a very useful probe in determining cosmolog- 
ical parameters. Because large and massive halos form very 
late, the high-mass tail of the mass function - the regime of 
cluster-scale masse s - is exponentially s ensitive to dark energy- 
related parameters dHaiman et alj|200lt) . Additionally, the red- 
shift evolution of the cluster mass function depends strongly on 
the cosmological parameters in a way that is complementary 
to other probes. The cluster mass function can also be used to 
measure the normalization of primordial density fluctuations, 
cr%, and search for h i nts of primo rdial non-Gaussianity (see, 
e.g jDalal et al.ll2008clOgurill2009l) . On cluster mass scales and 
smaller, the mass function, both directly and indirectly, plays an 
important role in halo models of galaxy formation and statistics, 



as applied to a wide range of redshifts and objects (predictions 
of bias, early galaxies, groups, quasars, spatial statistics of lu- 
minous galaxies, etc.). 

An important motivation for the precision determination of 
the mass function is the existence of several ongoing and up- 
coming surveys that aim to detect clusters via their optical, X- 
ray and Sunyaev-Zerdovich (SZ) effect signatures ( Rozo et all 



2010; Abbott et al. 2005; Vikhlinin et al. 2009; Kosowskv 2003; 
Staniszewski et al. 2009; Bartlett et al. 2008; Fang & Haiman 
2008). The number of detected clusters from the individual 
surveys will range from thousands to tens of thousands. To 
maximally extract cosmological information from these cluster 
surveys, the mass function must be specified to better than a 
fe w percent accurac y for a range of cosmolo gies. A s discussed 
bv lCunha & Evrardl j2010h and bv lWu et al.1 fcoid) , the current 
theoretical uncertainty in the determination of the mass func- 
tion can lead to a considerable degradation in the constraints on 
cosmological parameters. The investigations have also pointed 
to the usefulness of determining the mass function over a wide 
range of masses, extending down into the group scale. 

Because massive halos are very nonlinear and dynamically 
nontrivial objects, a fully satisfactory first principles approach 
to determine the structure and statistics of halos does not yet ex- 
ist. It follows, therefore, that our current theoretical understand- 
ing of the halo mass function is somewhat limited. (This situa- 
tion may be contrasted to that of nonlinear perturbation theory 
for the matter power spectrum, where independent of whether 
individual approaches fail or succeed, the actual problem is well 
defined conceptually and mathematically.) From the analytical 
standpoint, the only viable approach to the mass function is still 
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that ba sed on the (heuristic) Press-Schechter (PS) ex cursion set 
model dPress & Schechterlll974tlBond et al. II 19911) and its ex- 
tensions (see IZentner 1120071 for a review). Although this work 
has been valuable in suggesting functional forms and represen- 
tations for the mass function and in analyzing such effects as the 
scaling of finite-volume corrections, it has not independently 
yielded predictions for the mass function that are anywhere 
close to the accurac ies that are now requir ed. (For a recent criti- 
cal assessment, see Ro bertson et alj|2009l) . Moreover, it is hard 
to imagine how additional dynamics, gas physics, and feedback 
mechanisms can be modeled within such a framework. There- 
fore, it appears that a sufficiently accurate prediction for the 
mass function of halos can only be achieved using high reso- 
lution simulations, and modeling a range of physics tailored to 
specific applications. 

Numerical simulations have become a standard tool to de- 
termine the halo mass function over the last decade. Several 
groups have used suites of simulations to calibrate the halo 
mass fun ction over an increa s ingly wider range of masses and 
redshifts dJenkins et alj|200lt lEyrard et al.ll2002t IWhitej [2002; 
iReed et al. 1 120031: IWarren et al. Il2006t iHeitmann et al. [12006; 



rors, in the absence of rigorous testing of the universal ansatz 
they cannot be directly applied to cosmologies other than those 
considered specifically (and even in this case, the actual sys- 
tematic errors have often turned out to be larger than originally 
estimated). 

Motivated by these considerations, it is important to first es- 
tablish just how accurately various mass functions can be com- 
puted and what the systematic errors are in the most fundamen- 
tal situation - the gravity-only TY-body case. Once an accu- 
rate mass function for a particular ACDM case has been estab- 
lished, it is important to consider a range of observationally rel- 
evant redshifts and of cos mologies around that reference point 
(see, e.g., the discussion in lHeitmann et al.l2009l) . to understand 
and explore the range of applicability - and limitations - of the 
(almost) universal description described above. Therefore, the 
major aims of this paper are: (i) to carefully consider the sys- 
tematic effects due to numerical errors on the mass function 
and either avoid or correct for them; (ii) based on these results, 
establish an accurate prediction for the mass function of a ref- 
erence ACDM model at z = 0; and (iii) extend the investigation 
to a lar ger redshift range and provide an accurate prediction for 



Reed et al. 2007; Lukic et al. 2007; Tinker et al. 2008; Bovlan-Kolchimrelralgeneral vvCDM models (where the dark energy equation 



120091: ICrocce et al.ll2010h : see Uenkins et al.1 d2001l) for refer- 
ences to previous work. A key aspect of the calibration of 
the mass function is the use of In a~ l (M, z) as the central vari- 
able, instead of the halo mass, M. Here a 2 (M,z) is the vari- 
ance of the linear density field, extrapolated by linear theory 
to the redshift of interest, z, and smoothed by a spherical top- 
hat filter of radius R, which on average encloses a mass M 
(R = \3M/4Trpb(z)] 1 / i ). The associated scaled differential mass 
function f{cr,z;X) dJenkins et al.ll200ll) is 

f(a,z;X)=— (1) 

p b dln[a l {M 7 z)] 

where X labels the cosmological model and particular halo def- 
inition. The variable In a' 1 {M,z) appears naturally in the PS 
approach and extensions thereof, presenting a relatively simple 
form for f(a,z',X), in fact one with no dependence on cosmo- 
logical epoch and parameters. iJenkins et alj d2001l) found that 
for a certain fixed definition of halo, independent of cosmol- 
ogy, their simulation results covering redshifts from z = - 4, 
and across different cosmologies, could be well fitted by this 
"universal" form of the mass function at accuracies of order 
20%. Recent work has shown that mass function univers ality is 
apparently not valid beyond the 5 - 10% level dReed et alJ2007t 



Lukic et al.l2007l:lTinker et al.l2008tlCohn & Whitel2008ilCrocce et~ajf| 



2010tlCourtinetal.ll2010h 



Efforts to study this issue further quickly encounter a host of 
complications. Even independent of such significant physics is- 
sues as mass-observable mapping and baryonic effects, it turns 
out that the choice of halo definition and systematic errors in 
simulations can easily have as large an effect as that being in- 
vestigated. Thus, despite the major effort expended in numer- 
ical determination of the halo mass function, the present situa- 
tion cannot be considered to be fully satisfactory, as we discuss 
in Section [3] Among other sources of error, the effects of fi- 
nite force resolution and finite sampling error must be carefully 
dealt with in order to obtain a converged result. 

Beyond this point, there is a further cautionary note to keep in 
mind in terms of precision determination of the mass function: 
most simulation campaigns have focused on a single cosmol- 
ogy at a time. Therefore, even though results are often quoted 
in the universal form of Equation (JTJ with small statistical er- 



of state parameter, w, is constant in time, but w^-\). 

As a first step, the choice of halo definition has to be consid- 
ered. For the most part, numerical simulations use two different 
techniques to identify halos: friends-of-friends (FOF) or spher- 
ical overdensity (SO). In the FOF method, halos are found by a 
percolation technique where particles belong to the same halo if 
they are within a certain distance (the linking length b) of each 
other. The linking length is typically chosen between b = 0.15 
and b = 0.2, where b is defined with respect to the mean inter- 
particle spacing. The FOF definition of halos approximately 
traces isodensity contours and connects more directly to the 
simulated mass distribution; it is often used in cluster SZ stud- 
ies. However, the choice of linking length is an issue: too large 
a linking length can connect neighboring overdensities in a pos- 
sibly unrealistic manner. The SO method measures the mass in 
spherical shells around the center of the halo (which is usually 
determined from the potential minimum of the halo or from the 
most bound particle) until the density in the shells falls below 
a certain threshold which is given with respect to either criti- 
cal or background density. Typically, values for OT200 to wigoo 
(or higher in the case of clusters) are measured (with respect 
to p c ). The spherical overdensity method is particularly conve- 
nient for providing predictions for certain kinds of observations, 
lg. X-ray cluster masses where one is concerned primarily 
with studying the inner, virialized region of a halo. The major 
disadvantage of the SO method is the crudeness of the spherical 
approximation and that neighboring halos can overlap. 

Because isolated, relaxed halos are well-fit by the Navarro- 

SO and FOF 



Because isolated, relaxed halos are well-fit by 
Frenk- White (NFW) profile dNavarro et al.lfl997D. 



masses are strongly correlated dWhitel2001 



Tinker et al.l2008l) . 



In fact, if the halo concentrations are known, it has been shown 
that - in cosmological simulations - a one-to-one mappi ng for 
the tw o halo definitions exists at the 5% level of accuracy dLukic et al.1 
120091) . However, a fair fraction of halos in simulations are ir- 
regular. For currently favored cosmologies, 15-20% of b = 0.2 
FOF halos have irregular substructure or have two or mor e ma- 
jor halo components linked together dLukic et alj|2009l) . For 
such irregular halos, not only does the simple mapping between 
SO and FOF halos fail, it is not obvious just how to define an 
appropriate halo mass (lower b to what value, or correspond- 
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ingly, what choice of overdensity criterion to use?). 

In the absence of a compelling theoretical motivation, most 
numerical studies of the mass function have used FOF masses 



Table 1 
Parameters for the 38 cosmological models 



with linking length b - 0.2 following the convention set by Jenkins et 
d2001l) who noted that this definition led to a universa l form 
for the mass function (for a systematic investigation, see I White! 
l2002t lTinkeretal.1l2008l) . While noting its possible deficien- 
cies, we retain this convention here in order to better compare 
our results with other work. 

Our study of the FOF mass function uses a large suite of 
runs for a single reference A CDM cosmology (ver y close to 
the WMAP5 parameters from lKomatsu et al]|2009l) and a set 
of wCDM cosmological simulations - the "Coyote Universe" 
suite named after the supercomputer it was run on - that include 
a different ACDM model and a few others close to ACDM. 
The latter set of simulations represents a simple step beyond 
ACDM, where the dark energy equation of state parameter is 
treated in a purely phenomenological context. Allowing for 
dark energy evolution (as required by quintessence models, for 
example) opens up a large parameter space that near-future ob- 
servations are unlikely to be sensitive to. We, therefore, defer 
this extension to future work. 

In order to carefully control errors, we have followed the cri- 
teria for s tarting redshift, and mass and force resolution as pre- 
sented in iLukic et alj d2007l) . These criteria ensure that halos 
of a cer tain size and at certain redshifts can be resolved reli- 
ably. In lHeitmann et al.l J2009h similar criteria were laid out to 
obtain the matter power spectrum at 1 % accuracy out to scales 
k ~ 1 /iMpc" 1 . These criteria are also obeyed by the simulations 
used here. As discussed further in Section[3j the high-mass tail 
of the mass function is particularly susceptible to systematic 
errors in the determination of individual halo masses. These 
systematic errors can arise from the effects of finite force res- 
olution and we study and characterize these effects. Overall, 
the contribution of various errors in typical cosmological sim- 
ulations breaks d own basically as f ollows: (i) too low starting 
redshift, ~ 10% dLukic et alj|2007l) . (ii) halo sampling errors, 
~ 5%, (iii) transfer function approximations, ~ 5%, (iv) finite 
volume effects, ~ 1%, (v) force resolution effects, ~ 1%. Of 
these, (i) and (iii) are trivially avoidable, and the others can be 
controlled at least to the percent level. 

In this paper, after compensating for finite sampling and force 
resolution limitations, we present quantitative results for the 
mass function from the reference ACDM simulations and for 
a suite of wCDM cosmologies designed explicitly to bracket 
the curre ntly observationally rel evant range of cosmological pa- 
rameters dHeitmann et alj|2009l) . We provide a fitting formula 
describing our reference ACDM simulation data at the 2% ac- 
curacy level at the current epoch (we also compute the halo 
model prediction for the large scale halo bias from the mass 
function fit). We trace the evolution of the mass function be- 
tween redshifts z = 0-2, which represents up to a 10% break- 
ing of the universal description of the mass function across a 
representative range of masses. 

We then turn to investigating the variation of the mass func- 
tion as a function of cosmological parameters from the suite of 
wCDM simulations. This simulation suite, while it lacks some 
of the statistical power of the reference ACDM runs, provides 
a good test of the validity of the universal description of the 
mass function. At z = 0, we find that universality for differ- 
ent cosmologies holds to no better than at the 10% level, with 
clear systematic deviations (in both directions) from the quasi- 
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Note. — See text for more details and Heitmann et al. 12009) for the model 
selection procedure. To obtain good statistics over a wide range of halo masses, 
the reference ACDM case, model 0, was augmented by a set of additional runs 
(see Tablef2). The rightmost column shows the mass corresponding to \ja = 1.! 
for each cosmology at z = 0. 



universal form fitted to the reference ACDM runs. Thus, while 
not adequate for future precision studies (as is the case for all 
current mass function fits), our analytic form provides a good 
estimate for the mass function for wCDM cosmologies at the 
accuracy of currently available data. 

Our results demonstrate that the goal of determining the mass 
function at the percent level of accuracy will require a much 
more intensive program of simulations in the future, sampling 
both cosmological and physical modeling parameters (baryonic 
physics, feedback), along with w ell-controlled statis tical errors. 
As has been emphasized earlier dLukic et alj|2007l) . a possible 
solution is mass function emulation from a large, but finite set 



Mass Function Predictions Beyond ACDM 



of simulations, using techniques that have been s hown to be 
succe ssful in high-dimensional regression problems (lHabib et al. I 
l2007h . 

The paper is organized as follows. In Section[2]we describe 
the simulation suite used in this paper, encompassing 67 high- 
resolution simulations for 38 different cosmologies. Several 
overlapping-volume ACDM simulations are used to understand 
and control systematic errors. These errors and their ramifica- 
tions for the accuracy of the mass determination of halos, and 
how these translate to limiting the accuracy of the mass function 
itself, are discussed in Section [3] In Section [4] we present our 
results for the mass function for the reference ACDM model at 
different redshifts and provide a new fitting form for the mass 
function matching our simulations at the 2% level as well as 
the associated mass function-derived halo bias. We extend our 
discussion in Section [5] to the wider set of wCDM cosmolo- 
gies and investigate how well the mass function fit derived for 
the reference cosmology holds for this broader class of models. 
We conclude in Section [6] We discuss error control issues and 
provide relevant details in Appendix lAl 

2. SIMULATION SUITE 

Our simulation suite spans a wide range of observationally 
relevant wCDM cosmologies, as specified in Table Q] For each 
model we have results from a 1.3 Gpc box simulation, run with 
1024 3 particles, with masses of ps 10 10 M Q , exact values de- 
pending on the specific cosmology. We vary five cosmological 
parameters within the following boundaries: 

0.120 <u) m < 0.155, 

0.0215 <uj h < 0.0235, 

0.85 < n s <1.05, (2) 

-0.130 < w <-0.70, 

0.61 < rr 8 <0.90. 
The Hubble parameter h is fixed for models 1-37 by impos- 
ing the cosmic microwave background constraint, Ia = 7tdi s /r s = 
302.4, where d^ is the distance to the last scattering surface and 
r s is the sound horizo n. For a detailed descrip tion of the model 
selection process, see IHeitmann et al.l ( 20091). T he simulations 
are carried out with GADGET-2 dSpringelll2005l) . a tree-particle 
mesh (tree-PM) code. For a detailed discussion and compari- 
son of different A^-body methods use d for cosmological simu - 
lations, including GADGET-2, see, e.g. lHeitmann et al.ld2008al) . 
We use a 2048 3 PM grid and a (Gaussian) smoothing of 1.5 
grid cells. The force matching is set to six times the smoothing 
scale, the tree opening criterion being set to 0.5%. The soften- 
ing length is 50 kpc. 

For model 0, the reference ACDM cosmology, we have car- 
ried out additional simulations for different box sizes and mul- 
tiple realizations in order to cover a wide halo mass range be- 
tween 6- 10" M to 3-10 15 M with good statistics. Results 
from the overlapping volume boxes are also useful in under- 
standing systematic errors. Details about these ACDM simu- 
lations are given in Table [2] In addition to GADGET-2 we use 
a second tree-PM code for a subset of these simulations, de- 
scribed in I White! d2002l) . The algorithmic structure of this code 
is very similar to GADGET-2 a nd the code was also par t of the 
code comparison carried out in IHeitmann et al.l d2008al) . Aside 
from the main simulation runs, we also use a PM simulation 
with identical cosmological parameter settings as for the "G" 
run solely to study the impact of force resolution on individual 
halo masses. 



3. HALOS IN NUMERICAL SIMULATIONS 

An important question to understand before proceeding fur- 
ther is how uncertainties in individual halo masses - defined any 
way one chooses - translate into errors in the mass function it- 
self. To investigate this sensitivity we carry out some simple ex- 
periments. Introducing a Gaussian (or other symmetric) noise 
with a = 4% in individual halo mass measurements can make a 
small difference in the mass function at high masses (at the per- 
cent level) and should not be a concern. The picture changes 
substantially, however, if the halo masses are given small sys- 
tematic shifts. The error induced in the mass function can be- 
come much larger than the level of systematic error introduced 
in individual halo masses. We demonstrate this effect in Fig- 
ure [T| where systematic shifts in halo masses of 1%, 2%, and 
4% are studied. An increase of the masses by only 4% results in 
a difference in the mass function of 15% at high masses, which 
is quite significant. 

These results, arising from the exponential sensitivity of the 
mass function at high masses, demonstrate an important point: 
Given the level of uncertainties in halo masses due to numerical 
errors and the definition of the halo mass itself, it is not possible 
to derive a mass function prediction from simulations at sub- 
percent or percent accuracy without a consistent study of how 
individual halo masses vary as a function of simulation param- 
eters like force resolution, time step size, starting redshift, etc. 
Fortunately, it is already known that (b=0.2) FOF halo masses, 
over the mass range of interest, are relatively r obust to changes 
in sim ulation parameters; as demonstrated in IHeitmann et al. I 
d2005l) . individual halo masses as computed by six different 
codes with varying resolutions and time-ste pping schemes typ- 
ically agree to better than 2%. Additionally. iLukic etakl (120071) 
have provided criteria for running simulations so as to mini- 
mize systematic errors from a variety of possibilities. Our task 
is to ascertain whether error control can be further improved 
systematically. 

We focus here on three main possible systematic errors in 
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FIG. 1 . — Sensitivity of the mass function to systematic shifts in individual 
halo masses. Changes are shown relative to a baseline mass function, taken to 
be the fitting form of Tablef4] A small shift of 2% in the halo masses can lead 
to changes of up to 5-10% in the high mass tail of the mass function. 
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Table 2 
Specifications of the Simulation runs 










Box size 


Name 


n p 


m p 


„mm 
n h 


e 


Code 


£in 


Zout 


*- *runs 


ICs 


ACDM 


1000 Mpc 
1736 Mpc 
2778 Mpc 
178 Mpc 
1300 Mpc 


C 

B 
A 
GS 
G 


1500 3 
1200 3 
1024 3 
512 3 

1024 3 


l.MO lu M 
l.l-10 n M Q 
7.2-lO u M 
1.5-1O 9 M 
7.4-10 10 M Q 


400 

400 
400 
400 
400 


24kpc 
51kpc 
97kpc 
14kpc 
50kpc 


TreePM 

TreePM 

TreePM 

GADGET-2 

GADGET-2 


100/75 
100 

100 

211 
211 




0,1 

0,1 

0, 1,2 

0, 1,2 


2 
6 
10 
10 
2 


ZA/2LPT 
2LPT 
2LPT 

ZA 

ZA 


wCDM 


1300 Mpc 


Coyote 


1024 3 


varies 


400 


50kpc 


GADGET-2 


211 


0, 1,2 


37 


ZA 



Note. — Box size, mass and force resolution for the different runs; the upper section of the table describes the reference ACDM simulation suite (model in 
Table ff} while the lower section specifies the Coyote Universe runs (models 1 - 37 in Table ff). The total number of particles is denoted by n p , the particle mass by 
m p , n™" the nu mber of particles in the smallest halo kept, e the force resolut ion, and N mns the number o f realizations. For some simulations we used the Zel'dovich 
approximation IZel' dovich 197f|) (see also discussions in|Lukic et al. 2007 and Heitmann et al. 2010) to generate the initial condtions and 2LPT (Bouchet et al. I 
[J995lCrocce et aU2006D for others. 



the ma ss f unction (based on e xperience from iHeitmann et al. I 
d2005l) and iLukic et al] (120071) ): (i) small systematic errors in 
individual halo masses due to force resolution, focusing atten- 
tion on the high mass tail, (ii) the known systematic bias in in- 
dividual FOF masses as a function of the number of particles in 
the FOF halo dWarren et al. 1120061) . and (iii) systematic errors 
from missing long wavelength power due to the necessarily fi- 
nite box size. We also note the easily avoidable pitfall of using 
approximate fits for the transfer function rather than the numer- 
ical solution (see AppendixlAlfor details). As discussed below, 
and in more detail in AppendixlAl all of these effects can induce 
systematic errors in the mass function and need to be taken into 
account. (Figure[2]shows the effect of various corrections.) 

All numerical simulations are necessarily run with a finite 
force resolution, balancing the need for high spatial resolu- 
tion with suppression of noise/collisional artifacts. Below the 
chosen force softening length, the forces between the particles 
asymptote to zero. As a result, the halos that are formed tend 
to be "puffed out" in a simulation with coarse force resolution. 
The net effect is that - for heavy halos, with scale radii signif- 
icantly greater than the softening length - the mass of a given 
halo in a simulation with coarse resolution tends to be greater 
compared to a simulation with better resolution. This implies 
that the halo mass and hence the mass function determined from 
a simulation study would be higher compared to the "ideal" 
case of infinite force (an d mass) resolution . Although this ef- 
fect is known to be small (ILukic et al]|2007l) . it can certainly be 
significant at the percent level of error in the mass function. We 
also note that the effect of force resolution depends on the de- 
tails of how halos are found using the SO and FOF algorithms. 
Here we focus on the case of FOF halos; a disc ussion of the 
impac t of force resolution on SO halos is given in lTinker et al] 
d2008l) . In AppendixlAlwe provide a detailed analysis of the er- 
rors due to finite force resolution and how to correct for them. 
In our simulations, we find that finite force resolution effects 
can be accounted for by introducing a corrected rescaled mass 
M c via 

M c /M= 1.0-0.04(e/650kpc). (3) 

Here M is the uncorrected halo mass and e is the force reso- 
lution measured in kpc of the different runs as specified in Ta- 
ble |2] In our case, the biggest correction applied is w 0.6% for 
individual halo masses for the A runs (with a force resolution of 



97 kpc). This results in a systematic lowering of approximately 
2% in the high-mass tail of the mass function (primarily run A) 
as shown in Figure [AT6l in the Appendix. 

As stated earlier, we identify halos with a standard FOF algo- 
rithm, with a linking length, b = 0.2. Although halos with only 
a small number of particles (~ 20) can be reliably found with 
the FOF algorithm, accurate mass estimation requires keeping 
many more particles within individual halos. Aside from sim- 
ple considerations of particle shot noise, there is an inherent 
systematic error and scatter in the definition of an FOF halo 
mass with particle numbe r, even in the absence of all other lim- 
itations, as pointed out bv lWarren et al. I d2006l) . For ideal NFW 
halos (and for isolated relaxed ha los in simulations), this effect 
was studied bv lLukic et al] d2009l) and represents the best possi- 
ble scenario. To avoid problems with too few particles in halos, 
we restrict attention to halos with at least 400 particles. With 
this cutoff, the agreement in the mass function for the over- 
lap regions across the various boxes is within a few percent. 
After first app l ying t he FOF sampling correction suggested by 
IWarren et al. I (120061) . we find that a slightly modified correc- 
tion of the form «™ 1T = «/,(l -«^ 0,65 ) brings the results from the 
nested boxes in good agreement, as shown in the Appendix in 
Figure IA17I We stress that this adjustment is purely empir- 
ical, targete d at matching halo masses in overlapping boxes. 
The study in ILukic et all (120091) shows that for NFW halos, the 
correction depends on «/, as well as on the halo concentration. 
Thus, in principle, one may expect the FOF sampling correction 
to be more complex than a simple compensation based on «/,; 
we leave a more detailed analysis for future work. The use of 
the restriction «/, > 400 appears, however, to make the correc- 
tion predominantly dependent on «/, alone. The net correction 
in halo mass due to finite force and mass resolution, for our 
simulation suite, can then be written as 

M c /M= [1.0-0.04(e/650kpc)](l-n^ - 65 ). (4) 

Last, we consider systematic errors due to the finite volume 
of simulations. There are three sorts of effects of this type. 
The first is simply that the number of halos at high masses will 
be poorly sampled, and the mass function in this region will 
have large statistical error bars due to shot noise (Poisson fluc- 
tuations). This is purely a question of having sufficient total 
simulation volume. The second effect is the fact that missing 
large-scale modes, with k < 2-k/L (L is the box size in lin- 
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FIG. 2. — Uncorrected (upper panel) and corrected (lower panel) results for the ratio of the mass function to our best fit model at z = 0. Most of the difference is 
due to the correction for the finite particle sampling of the halos. A detailed description of the error bars is given in the appendix. 



ear dimension), lead to a suppression of structure formation, 
and hence of the mass function, as reduced power is available 
for transfer from linear to nonlinear scales as evolution pro- 
ceeds. Mass functions measured from simulations must there- 
fore account for the infrared cutoff in the variance of matter 
fluctuations a(M). The extended Press-Schechter approach has 
been found to wo rk well in compensating for this effect (see 
iLukic et alj|2007l and references therein), and we follow it here. 
For our simulations, this volume correction is relevant only for 
the small box (GS) set of simulations, where L = 178 Mpc, and 
affects only the low mass halos. Applying the EPS method to 
the smallest box shifts the masses by 2%. 

The third finite volume effect is related to the (effective) num- 
ber of independent realizations, i.e., the sample variance. Be- 
cause halos are biased tracers of the density field, the mass func- 
tion in a given, sufficiently large, target volume is sensitive to 
the local mean density as set by the scale-independent bias. In 
a large-volume simulation, local sub-volumes will have fluctu- 
ating average densities and these will add another component 
to the mass function variance, aside from shot noise. There 
are two ways to address this: (i) because of the high covari- 
ance between the small number of low k modes in a simulation 
and its high k modes, run either very large boxes where the 



relevant low k modes (in terms of power transfer) are sampled 
sufficiently well, or (ii) run a sufficient number of statistically 
independent large-volume boxes. 

In either case, one can successfully estimate the variance us- 
ing a simple halo model prescription and linear theory for den- 
sity fluctuations dHu & Kravtsov 1 120031) : details are given in 
Appendix |A.3l An important point to note is that all simulation 
boxes necessarily have a definite infrared cutoff set by the box 
size. This can in fact be tuned to control the mass function vari- 
ance: smaller boxes have smaller variance because they have 
less low-frequency power. Of course, one cannot make the box 
too small because then the mass function will be biased low as 
in the EPS discussion above. 

If one runs only one box, then resampling techniques such 
as the jackknife may have to be used to estimate the associ - 
ated errors (see, e.g., iTinker et al.l l2008t ICrocce et al.1 120101) . 
These techniques can be susceptible to misestimation of errors 
for smaller boxes, due to the assumed independence of subvol- 
umes of the simulation box. To estimate our errors, we prefer 
to use independent realizations for a given simulation volume. 
This has the advantage that the individual modes of fluctuations 
are truly independent across realizations and hence averaging 
over them represents the true variance for each set of runs. We 
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note that for large b oxes ^2 Gpc, jacknife sampling works well 
dTinker et al.ll2.Q08b . However, since we have 10 independent 
2 Gpc box runs (A), we can use the different realizations to es- 
timate the errors. (Additionally, the increased mass resolution 
in each box makes one less susceptible to the finite sampling 
problem in determining FOF halo masses.) The sample vari- 
ance errors were computed by taking the ratio of the variance 
for each of the runs divided by the number of realizations of 
each run; details are discussed in Appendix |A. 31 As stated al- 
ready, with this method one also has the extra flexibility of tun- 
ing the low frequency cut-off to optimize the mass variance (if 
the only simulation goal is to produce an accurate mass func- 
tion). 

Our results from this section are summarized in Figure [2] 
The upper panel shows the raw simulation results and the lower 
panel includes our three corrections for force resolution, finite 
sampling, and finite volume. In order to show the effects more 
clearly, we display the ratio of the mass function to our best 
fit to the data, as derived in the next section. Once the simu- 
lation parameters for initial redshift, force resolution, and box 
size are chosen in a suitable way for studying the mass func- 
tion, the major correction required is the one due to finite mass 
resolution. 

4. REFERENCE ACDM MODEL 
4. 1 . Mass Function at the Present Epoch 

After having analyzed possible systematic errors and cor- 
rected for them in all our simulations, we now investigate the 
mass function for the reference ACDM cosmology specified in 
Table Q] (model 0) at z = 0. The effective simulation volume, 
combining all of our runs, is approximately 250 Gpc 3 . The 
simulations provide coverage of halo masses ranging from that 
relevant for individual bright galaxies all the way to clusters, 
with good statistics. Although the total volume is dominated by 
the set of runs with box size of 2778 Mpc, the other runs are 
very helpful in checking for and investigating systematic errors 
as seen above, and in extending the simulation reach to lower 
halo masses. 

As previously discussed, a convenient form to express th e 
scaled differential mass function f(cr,z) is dJenkins et al.l200lh : 
dp/p b _ M dn(M,z) 
J(<T,Z> dlna- 1 p b (z)dln[a- l (M,z)]' 
Here n(M,z) is the number density of halos with mass M, ph(z) 
is the background density at redshift z, and a(M,z) is the vari- 
ance of the linear matter power spectrum P(k) over a length R, 

D\z) 



a z {M,z). 



2tt 2 



k 2 P{k)W 2 (kR{M))dk, 



(6) 



when smoothed on the scale R(M) = QM/A-irpt,) 1 ' 3 with the top- 
hat filter W(x) = 3[sin(x)-xcos(x)]/x i . We write a{M) = a for 
brevity in the following. The redshift dependence is encapsu- 
lated in the growth factor D(z) which is normalized in such a 
way that D(0) = 1 . As mentioned earlier, the advantage of this 
definition of the mass function is that to a reasonable accuracy 
it does not explicitly depend on redshift, power spectrum, or 
cosmology; all of these are encapsulated in a(M,z). 

A popular nu merical fit for th e diffe rential mass function 
f(a) is given in ISheth & Tormenl d 19991) (ST hereafter). The 
expression for the ST mass function is 



fsr(o-) 



-A\ — exp 

V 7T 



aSl 
'la 1 



1 + 



id 2 



-, (7) 



where A, a, and p are three parameters tuned to simulation re- 
sults with a = 0.707 (a = 0.75 is proposed as a better estimate in 
ISheth & Tormenll2002l ). p = 0.3, and A = 0.3222. The parame- 
ter A is fixed by the normalization condition that all dark matter 
particles reside in halos, i.e. 



dlnaf(cr) = 1. 



(8) 



We note that, as a practical matter, the lower halo mass cut-off 
in numerical simulations is too large to test this particular as- 
sumption. So, in principle, one could leave this constant as a 
free variable in the fitting process. With the normalization con- 
dition fullfilled, the ST mass function has two free parameters, 
a and p. 5 C is the density threshold for spherical collapse. In an 
Einstein-de Sitter cosmology, S c = 1.686, independent of red- 
shift. For n,„ y 1, the value for 5 C shows insignificant depen- 
dence on cosmology dLacev & Coldll993l) . We checked that 
including the cosmology dependence in 6 C does not explain the 
redshift evolution of /(er) seen in our simulations. In the fol- 
lowing, we therefore keep S c = 1 .686 as a fixed value. 

In order to obtain a fit for f(a,z), we need to compare Equa- 
tion ([5]) with the binned mass function obtained from the sim- 
ulations. We measure the number density of halos in a bin of 
size AlnM with mass limits [M\,M2], as 

f I din a 
"bin = Pb / ———fdau(M,z)dlnM, 



AlnM 



Md\nM J 



(9) 



In the limit that AlnM —} 0, we can write Equation (O as 



«bin = Pb- 



1 dlna(M hm ,z) 



M bin d\nM hin 
and hence /data(Mbi n ,z) can be written as 

"binM bin (d In <r(M b 



/data(M bin ,z)AlnM bil 



/data(M bin ,z): 



1.2 



p b AlnM hil1 \ <fln M b 



i,z) 



(10) 



-*-bin= 242 
T -.-bin= 63 
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FIG . 3 . — Ratio of the mass function data at z = with respect to our z = fit 
for different choices of the number of bins. For halos of mass < 10 14 MQ,four 
bins per decade in mass are sufficient for the mass function data to converge. 
For halos of mass > 10 14 Mq, the result converges for a choice of 15 bins per 
decade. 
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Fig. 4. — Ratio of various mass function fits derived in previous studies with 
respect to the results of this paper at z = 0. The binned numerical data are the 
points with error bars; the ratios are ta ken with respect to the analytic fit to the 
numerical data specified by Equation H2\ . Because the fits are based on runs 
with different cosmologies, exact correspondence cannot be expected (since 
the mass function is not universal). 



where M\,[ n = ^Mi/N^v, and the summation is over all the ha- 
los in a bin. When combining all the boxes from the various 
simulations, we vary AlnM and make it small enough to en- 
sure that Equation (TTOb holds within the accuracy of the data. 
The results are shown in Figure [3] We find that four bins per 
decade equally divided in ln-space between 10 11 — 10 14 M and 
15 bins per decade between 10 14 — 3 • 10 15 M Q are sufficient 
for the mass function to converge within the measurement er- 
rors for z = and z = 1. This corresponds to a bin size of 
AlnM = 0.25 between 10 11 - 10 14 M© and AlnM = 0.06 be- 
tween 10 14 - 10 16 M . For z = 2, four bins per decade in mass 
(AlnM = 0.25) are sufficient for the measurements to converge 
within the accuracy of the data. 

Figure|4]compares different mass function expressions given 
previously as compared to our new simulation results at z = 0. 
The expressions for different fits are given in Table [3] Some 
of these expressions, derived from earlier simulations, are sig- 
nificantly discrepant - especially at high masses - at the 20% 
level and higher. Results from more recent simulations are in 
much better agreement, partly because the cosmologies run are 
much closer. We note that exact correspondence cannot be ex- 
pected because the mass function is not universal; we defer a 
discussion of this point to Section[5] 

In order to obtain a fitting form to our results, we begin with 
the ST fit as the starting point, although, as shown in Figure [4] 
the ST mass function deviates from the simulation results by as 
much as 40% at the high mass end. As a first step to improve the 
accuracy of the ST mass function, we drop the normalization 
requirement and refit all three parameters A, a, and p to the 
numerical data. With this approach, the simulation data can 
be fit to an accuracy of 10-15%, signific antly worse than th e 
statistical errors of our data set (see also iManera et alj|2010t) . 



Table 4 

Mass function fitting formula for the reference 

ACDM MODEL (MASS RANGE: 6 X 10" Mq- 3 x 10 15 M ; 

REDSHIFT RANGE: z =0-2) 



r°V,z)=Aj^exp 



a5 c 

~2o- 



1 + 



(*)' ( 



ScVa 



Redshift Evolution 



7 0.333 ~ 

n ~ (1+7)0.11 ; " ~ 



0.788 ~ _ 0.807 - _ 1.795 
(1+,-) " 1 , P ~ (1 +; )0.0 i H ~ (l+;)Q.o 



The remaining inadequacy o f the ST fit can be a ddressed in 
different ways. For example. IWarren et al. I d2006l) introduced 
a fourth parameter into the ST functional form and refitted the 
other three parameters to their simulation data, finding a best fit 
mass function: 



f w (<j)=A v 



1 



+ c I exp 



(11) 



with A w = 0.7234, b = 1.625, c = 0.2538, and d = 1.1982. As 
shown in Figure [4] this particular fit also severely underesti- 
mates the mass function at high masses, by up to ~ 30%. While 
adequate as a fitting form over a finite range, Equation ( fTTT i 
diverges when the normalization condition is imposed [Equa- 
tion dHJ]. To avoid this, we present a new fitting function for 
/(ct). This is the simplest ST modification that does not diverge 
but adds one extra parameter, qo (for qo = 1 we recover the ST 
mass function): 



/ mo V,z = 0)=A 



exp 



Q(A 2 

"2a 2 



1 + 



a Q 5- 



i>" 



Sa- 



il!) 

We use a x 2 technique to determine the best fit f(a) that matches 
the mass function data obtained by combining all of the ACDM 
runs. That is, we minimize 



X 



=E 



(A/(<T) data ) 2 



(13) 



where /(er) mod , /(er)data and A/(cr)d a ta are given by Equations (Tl2b . 
( [Tol l, and dA4l ) respectively. 

Minimizing x 2 gives the best fit parameter values: Aq = 0.333, 
do = 0.788, po = 0.807, and qo = 1.795 with a x 2 P er degree of 
freedom of 1.15. The subscript "0" indicates that the best fit 
values are specified at z = 0. The results are summarized in 
Table [4] As mentioned above, this expression does not diverge 
when the normalization condition is imposed, however, the best 
fit does not lead to a normalization of unity. As shown in Fig- 
ure [5] this modified expression agrees with the simulation data 
to better than 2% accuracy at z = 0. As further discussed in Sec- 
tion [4J2] a simple redshift dependence has to be introduced into 
the fitting function to obtain agreement at the same accuracy 
level at higher redshifts. 

4.2. Redshift Evolution and Universality 
The z = mass function fit of Section[4~T|has a default univer- 



sal form. However, as remarked previously, it is known that the 
mass function deviates from universality - as a function of red- 
shift - for ACDM cosmologies. Recent results include those 
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Table 3 
Mass function fitting formulae derived in previous studies 



Reference 



Fitting function f(er) 



Mass Range 



Redshift range 



Sheth & Tormen (2002) f ST (a) = 0.3222- '^™> 



^exp 



0.75(5; 



~la^~ 



1 + 



^oTsSjJ 



0.3 



Jenkins et al. (2001) 
Warren et al. (2006) 



Reed et al. (2007 



0.3222 



0.315exp [-liner" 1 +0.61 1 3 - 8 ] 
0.7234 (a" 1625 + 0.2538) exp [-1^2] 

(otj) ' ' + 0.6G l (a) + 0AG 2 (cr) 

0.03 



2(0.707) 



1 + 



:-exp 



Manera et al. (2010) 
Crocceetal.(2010) 



0.764(5,; 

2a 1 



(n eff +i)\S c laf 



fsT(a) = 0.3222, f exp 



~2(T 2 



1 + 



(*)' 



A(z)[<J- a ^ + b(z)]exp[- c -£] 



Unspecified 
-1.2 < In a- 1 > 1.05 



(10 10 - 10 15 ) h^Mc 



-0.5 < In cr 1 > 1.2 



(3.3 x 10 13 -3.3x 10 15 )h-'M G 



(lO^-lO'^h-'Mg 



Unspecified 

z=0-5 
z=0 

z=0-30 



z=0-0.5 
z=0-l 



Note. — Various fits from previous studies shown in Figure f4] and [6] for friends-of-friends halos of linking length b = 0.2 are listed. For Manera et al. (2010), 



the parameter values are (a, p)= (0.709, 0.248) at z=0 and (0.724, 0.241) at z=0.5. For lCrocce et all 42O10D . the parameter values are A(z) = 0.58(1 + z)~ 



,a(z)- 



1.37(1 +z)-°- 15 ,fe( z ) = 0.3(1 +zT 



,c(z) = 1.036(1 +zf 



For lReed et al] UMfo , G\(o) = exp 



2(0.6) 



-0.4)- 

~7— 



and G2(cr) = exp 



(ln<T~'-0 75)- 
2(0.2) 2 



of iTinkeret alj d2008l) who found the SO halo mass function 
to evolve by 20-30% from redshift z=0 to 2.5. As shown in 
Figure [5] this deviation can be as much as 10% between red- 
shifts z = — 2 for F OF halos, in agreement with the results of 
ICrocce et alj d2010l) . In this section we extend our fitting func- 
tion to include the redshift evolution of the mass function. We 
parameterize the possible redshift evolution of each parameter 
via a simple power-law form 

A=A /(l+z) ai , 

s = a /(i+z) a2 , 

p = p /(l+z) a \ 

q = q Q /(l + z) a \ (14) 

In order to ensure that the expression for the redshift evolution 
reproduces the mass function at any intermediate redshift when 
interpolated or even extrapolated, we fit two redshift outputs at 
a time. Thus we have three values for each parameter. The fi- 
nal set of parameters is the average of the three values obtained 
using redshift outputs in pairs. Figure [5] shows that the power 
law model of Equations ( TBI is able to capture the redshift evo- 
lution with an accuracy of better than 3% within the range of 
0.6 < 1/cr < 2.4. We find that only two of the four parame- 
ters of Equations ( TBI show any redshift evolution. The best 
fit values for the parameters a, describing the redshift evolu- 
tion are a\ =0.11, a 2 = 0.01, aj = 0.0, and cxa, = 0.0. We note 
that the non-universal redshift evolution is suppressed at higher 
redshifts as the effect of the cosmological constant is reduced 
and matter-domination takes over. To recap, our analytic best- 
fit to the mass function data uses one extra shape parameter 
compared to ST to match the z = data, and then introduces 
a simple z-dependence (two more parameters) to capture non- 
universal behavior. 

We also explore the option of allowing for a redshift depen- 
dence, S c (z), as determined by the spherical collapse model. 



For the ACDM cosmology adopted here, the spherical collapse 
calculation yields S c = 1.674, 1.684, and 1.686 respectively for 
z = 0, 1, and 2. As expected, at z = 2, S c approaches the value 
expected for the £!„, = 1 cosmology. Allowing for this redshift 
dependence does not mitigate the redshift dependence of the 
other fitting parameters; including S c (z) changes the value of 
a to 0.799/(1 +z) 0024 (in fact increasing the redshift depen- 
dence) while the other parameters remain unchanged. There- 
fore, adding S c (z) does not help explain the z-dependence seen 
in our simulations, and we do not include it in our fit. 

High-statistics studies of the evolution of the FOF mass func- 
tion have been carried out previously. In an investigation focus- 
ing m ainly at high redsh ifts, to explain the violation of univer- 
sality, iReed et all J2007I) proposed an effective spectral slope 
« e ff set by the halo radius, parameterized as 

diner -1 „ ,„_ 

aflnM 
This new effective slope induces a redshift dependence in the 
m ass function . How ever, as shown in Figure [6] the analytic fit 
of lReed et al. J2007I) i s not i n good agreement with our results 
(see also iBagla et all d2009l) which studies the dependence of 
the mass function universality on the initial matter power spec- 
trum). This discrepancy indicates that high redshift evolution 
of the mass function is slo wer compared to that at lower red- 
shifts. ICrocce et alj d2010l) also use a simple power-law form 
to fit for redshift evolutio n. Ou r reference model and the pa- 
rameters of ICrocce et al.1 d2010h are very close, consequently, 
their results are significantly closer to ours, except at very high 
masses, where their fit appears to deviate. Here, discrepancies 
may result from the use of an approximate transfer function Q 
and a small systematic offset in their fitting procedure at high 
masses (Cf. their Figure 8). Taking this offset into account 

1 M. Crocce, private communication 
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Fig. 6. — Redshift dependent mass function fits as introduced bv lReed et al.l 
12007) and Crocce et al. 1 2010) compared with the numerica l data of thi s work . 
Aside from disagreement in the overall shape, the results of Reed et al. 12007) 
underestimate the amount of evolution indicating that high redshift evolution 
of the mas s function is s l ower compared to that at lower redshift. The agree- 
ment with ICrocce et al.l 120101) is better (at the 4-5% level), except for the 
runaway at high masses (see discussion in Section l4~2l . 



we find that the actual numerical results are in agreement at 
th e few percent leve l. Th e expressions for the fitting functions 
of lReedetail d2007l) and lCrocceet al.ld2010l) are given in Ta- 
ble[3] Figure [7] shows the abundance dn/dlnM as measured in 
our simulation along with the analytic fits. Figure [8] summa- 
rizes the results from this section, showing the mass function at 
different redshifts and our best fit results. 

4.3. Mass function-derived large-scale halo bias 

The evolution o f the spatial distributi on of halos has been 
studied in det ail inlCole & Kaiset (1 19891) and subsequently in 
iMo & White! (1 19961) and lShefh & Tormenl (Il999h . These stud- 
ies assume that dark matter halos are biased tracers of the un- 



derlying dark matter distribution. The halo bias in general is 
stoch astic and a nonlinear function of the underlyin g density 
field (ISchulz & White! [20061: ISeliak & Warred 12001 . In ad- 
dition, it also dep ends on the asse mbly history of the halo s 
JDalaletal.1120081) . As discussed in ISheth & Tormenl dl999l) . 
within the halo model, the large scale deterministic bias can 
be derived knowing the shape and evolution of the mass func- 
tion. It is known that the peak-background sp lit model pre- 
diction for the bias using the Sheth-Tormen or IWarren et al. I 
d2006l) mass function is in disagreement with direct numeri- 
cal calculations of this quantity obtained from ratios of correla- 
tion functions or power spectra in simulations (see, e.g- jLukicl 
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l2008tlPadmanabhan & Whitd2009l and lTinkeret alj2010h . Given 
a more accurate mass function, it is easy to check if the halo 
model approach now p roduces a better answer for t he halo bias. 
To test this - following lSheth & Tormenld 19991) and lCole & Kaiser! 
dl989l) - we now obtain the expression for the large scale bias 
using our expression for the mass function. 

The expression for bias can be written in terms of the condi- 
tional and unconditional mass functions as 

N(m,zi\M,V,Zo) 



b h 



n(m,zi)V 



1 



(16) 



"i/O^i) 

where N(m, z\ \M, V, zo) is the average number of halos of mass 
m which collapsed at z\ and in a cell of volume V which con- 
tains the mass M at zo and v = 5 2 /cr. 

In the large scale limit, the 'peak-background split' formal- 
ism dSheth & Tormenll999l:IC~ole & Kaiserl 19891) prediction for 
the halo bias for mass m at redshift z can be derived as fol- 
lows: Define the peak height v\ relative to the background v$ 
(&\ — 5o) 2 /(g 2 —Oq). Keeping the leading order terms in 



as v 



10' 



the expression gives v 2 = v 2 (\ — 28q/8\). One can then Taylor 
expand Equation dT6l ) and use Equation (01 to obtain the bias 
predicted by the ST mass function in Lagrangian space. Con- 
verting to bias in Eulerian space, the result is 



b S T = 1 + ■ 



1 



2p/S c 



(17) 



S c \+{av)P 
Similarly, using the mass function fit based on our simulation 
results, Equation ([T2I 1. the large scale bias can be expressed as 



''mod " 



1 + 



av-q 2p/S c 



(18) 



S c \ + {av)i>' 
Figure [9] shows the predicted bias compared to the ST result. 
Due to the difference in the two mass function fits, we find 
a corresponding change in the large scale bias of 10- 15%. 



Note that this difference is maximal for halos in the mass range 
10 13 - 10 14 M . Unfortunately, the fact that the (nominally) 
improved prediction for the bias lies consistently below the ST 
prediction, makes it deviate even further from numerical calcu- 
lations for the large-scale bias as a function of mass. This be- 
havior is in excellen t accord with the ana lysis in lLukicI Q2008) 
performed using the I Warren et al. I ( 20061) mass function fit and 
simulation data from iLukic et al.l d2007l) . Therefore, we con- 
clude that the simple halo model result for the halo bias does 
not converge correctly as one essential ingredient - the mass 
function accuracy - is systematically imp roved. Our conclu- 
sion is consistent with other recent studies: iTinker et al.l d2010l) 
have found that SO halos are systematically biase d 10% higher 
compa red to the Sheth-Tormen bias prediction. iManera et al.l 
d2010l) also conclude that a prediction based on a simple peak- 
background ansatz is inadequate to find agreement with their 
simulation results. 

5. MASS FUNCTION FOR wCDM COSMOLOGIES 

Based on our results in Section|4] which include a mass func- 
tion fit at 2% accuracy for the reference ACDM cosmology, we 
now investigate how well this fit works as cosmological param- 
eters are varied. Including the dark energy equation of state pa- 
rameter, w , as a phenomenological constant, we consider a suite 
of wCDM models: 37 simulations that span a wide range of cos- 
mological parameter values. The values of the parameters for 
each simulation are given in TableQ~]with the ranges being spec- 
ified in Equation (f2|. Although each simulation is individually 
large, the volume of each run is much smaller (2.2 Gpc 3 ) than 
for the combined reference ACDM runs (250 Gpc 3 ). Therefore, 
the statistics of the halos, especially at high masses, is not as 
well-determined in this case. Consequently, we do not attempt 
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Fig. 10. — Test of universality of the mass function for wCDM cosmologies (the G runs) at z=0. The range of the cosmological parameters covered by the 
simulations is given in Table [T] The ratio is taken with respect to the best fit mass function for the ACDM case at z=0. The lines show the ratio of the fit if S c is 
cosmology dependent in the ACDM fit. 



an as careful statistical fit for the mass function for the vvCDM 
models, although we do use all the relevant error controls. For 
each run we compute the Poisson error and the sample variance 
error using Equation A5. 

Interestingly, it turns out that the modified universal fit al- 
ready derived for the ACDM case, including the redshift evo- 
lution, also works in a relatively unbiased way for the suite of 
wCDM cosmologies at the 5- 10% level. However, clear sys- 
tematic violations of universality are observed from this suite 
of runs as cosmological parameters are varied (Figures [TOj \TT\ 
and[T2l. The figures show the ratio of the wCDM mass func- 
tions at z = 0, 1, and 2 obtained from our simulation suite with 
respect to the reference ACDM mass function. We find that 
the FOF mass function shows systematic ~5-10% variations (in 
both directions) with respect to the ACDM reference (at higher 
masses, the statistics is not good enough to make very precise 
statements), although universality is preserved at the 2-3% 
level upto 1/cr ~ 1.2. Note that for some choices of cosmo- 
logical parameters the difference is as large as 10% especially 
at the high mass end at z = 0. 

We also investigate the cosmology dependence of S c as pre- 
dicted by spherical collapse. Details of the calculation of S c as 
a function of cosmology are given in Appendix [B] (See also 



iPace et aUuOlOl for a detailed study.) The lines in Figures 10, 
11, and 12 represent the ratio of our fit at a particular redshift 
when S" CDM is used in our expression in Table 4 to the refer- 
ence ACDM fit at that redshift. Adding 5"' CDM agrees qualita- 
tively with the wCDM runs, however the difference predicted 
by the spherical collapse model is too small to be in quantita- 
tive agreement with our results. Note that S"' CDM depends only 
on two parameters for a flat wCDM cosmology, namely tt m and 
w. However, universality can also break down as other parame- 
ters change, and some of these - erg, h, and n s , - are also varied 
here. There is currently no theoretical framework to understand 
the breaking of universality with these parameters. At higher 
redshift, wCDM cosmologies are expected to converge to the 
fl m =l cosmology. This is roughly the trend seen in Fig. 1 1 and 
12 for the case of z = 1 and z = 2. 

Figure [13] summarizes our findings in this section showing 
the best fit mass function /(er) at three redshifts and the simula- 
tion results from the 37 cosmologies. Overall, all cosmologies 
are described reasonably well by our new fit. Points to be noted 
include: (i) An approximately universal fit describes the wCDM 
mass function at 10% accuracy over an observationally useful 
range of cosmological parameters, and (ii) universality of the 
mass function is systematically broken at this level of accuracy 
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Fig. 1 1 . — Ratio of the mass function for wCDM cosmologies to the best fit mass function for the ACDM case at z=l. The lines show the ratio of the fit if S c is 
cosmology dependent in the ACDM fit. 



within both wCDM (for example, models 22, 25, and 35) and 
ACDM cosmologies (models 31 and 34). 

Given the current observational state of the art, a wCDM 
mass function determined at the 5- 10% level over the range 
of cosmologies studied here appears adequate for data analysis. 
However, this is not likely to be the case for the next generation 
of observations. To extend our work further would require a se- 
ries of large simulations with their mass function results inter- 
polated along the same lines as already a chieved for the power 
spectrum to scales of order k ~ 1 /zMpc -1 (Lawrence et al.l2010l) . 

6. DISCUSSION 

In this work, we study the mass function of dark matter halos 
over a wide range of masses (6 • 10 11 - 3 ■ 10 15 M Q for a refer- 
ence ACDM model) and a redshift range of z = 0-2 for a large 
range of wCDM cosmologies. A friends-of-friends algorithm 
with a linking length of b = 0.2 is used for halo identification. 
The primary aims are to control numerical errors and gain suf- 
ficient statistical power for cosmological parameter estimation 
and testing of the universality of the mass function. For a refer- 
ence ACDM model, we achieve a 2% error in determining the 
mass function. At this level of numerical control, deviations 
from universality (in both cosmological parameters and redshift 
evolution) can be studied systematically. Our level of numeri- 



cal control is approaching the A^-body baseline level - neces- 
sary but by no means sufficient - required by next-generation 
cosmological surveys. The quest for high accuracy in the 'N- 
body' mass function has a natural stopping point at the percent 
level simply because at this point many oth er physical processe s 
become important (e.g., baryonic effects, Stanek et al. 2009). 
Moreover, the connections to observations need to be directly 
modeled and end up adding their own significant contribution 
to the overall error. 

We use a large number of high resolution simulations for 
studying the mass function in a single ACDM cosmology. These 
simulations are used to establish the error control methodology 
in order to obtain an accurate mass function. Error sources sys- 
tematically studied here include effects of finite force resolu- 
tion, FOF particle sampling bias, and systematics induced by 
finite- volume effects. 

Using the more accurate mass function fit for this reference 
cosmology, we also rederive the large scale halo bias using the 
'peak-backgound split' approach. Compared to previous re- 
sults obtained with the ST fit, the halo bias changes by 10-15% 
between z = - 2, and is systematically lower. Unfortunately, 
instead of improving the agreement with direct numerical mea- 
surements of halo bias (computed by taking ratios of correlation 
functions), the use of an improved mass function only increases 
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FIG. 1 2. — Ratio of the mass function for wCDM cosmologies to the best fit mass function for the ACDM case at z=2. The lines show the ratio of the fit if S c is 
cosmology dependent in the ACDM fit. 



the discrepancy. This points to an essential difficulty with the 
halo m odeling approach . We will return to this problem else- 
where dLukic et al.ll2010l) . 

After studying the mass function in detail for the reference 
ACDM cosmology, we extend the range of cosmological pa- 
rameters by considering a suite of wCDM simulations. The 
range of parameters is set by the current constraints on cosmo- 
logical parameters. The simulation parameters for the runs are 
given in Tables Q] and [2] We note that the reference mass func- 
tion fit provides a good description of the wCDM results at an 
accuracy of ~ 10%, however, with systematic deviations that 
point to clear violations of the universality of the mass func- 
tion, not only in wCDM parameter space, but also within the 
set of ACDM models (and models very close to ACDM). 

The breaking of univ ersality as a function of different ACDM 
parameters is s tudied in lJenkins et all (1200 ll) : lTinkeret al.ld2008l) : 
IWarren et al. I d2006l) . These studies have shown that universal- 
ity b reaks down at the 20 % level when a% and J7 m are varied. 
Also lCourtin et al.l ( 120101) have studied how universality breaks 
down for quintessence cosmology and found similar variation 
of universality with cosmology as reported here. The univer- 
sality of the mass f unction in the presence of early dark energy 
has been studied in lGrossi & Springe! d2009l) . We find that uni- 
versality in the mass function holds at the 10% level as a func- 



tion of cosmological parameters (the wCDM suite) and redshift 
(wCDM and the reference ACDM model). 

Given the break ing of universality discussed here (see also 
iTinker et al.ll2008l) . it is not clear what the utility of fitting for- 
mulae for the mass function might be for observations that actu- 
ally do require theoretical predictions at the percent level of ac- 
curacy. An extensive future simulation campaign will be needed 
to properly come to grips with this problem. Nevertheless, to 
provide a compact description of our results, we derive a fit- 
ting formula that agrees with our reference ACDM simulation 
results at 2% accuracy over the redshift range of z = 0-2 and 
the mass range of 6 ■ 10 1 ' - 3 • 10 15 M . The fitting function is 
a simple modification of the Sheth-Tormen form with one extra 
shape parameter to improve the fit at z = and two extra evolu- 
tion parameters. This form does not lead to divergences if the 
normalization condition that all mass resides in dark matter ha- 
los is imposed. It holds at the 10% level of accuracy for a broad 
class of wCDM cosmologies with the redshift evolution taken 
into account. 
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APPENDIX 

SYSTEMATIC ERRORS IN MASS FUNCTION MEASUREMENTS 
FROM N-BODY SIMULATIONS 

Initial Conditions 



1.15 



1.1 



C 



J? 1.05 



One issue regarding initial conditions not considered in lLukic et all 
d2007l) is the accuracy of the transfer function used for gener- 
ating the linear power spectrum. Usually computed via Boltz- 
mann codes such as CMBFAST dSeliak & Zaldarriagal[l996l) 
or the related code CAMB3 t he transf er function is in princi- 
ple known to 0.1% accuracy dSeliak et al. 112003b . quite suffi- 
cient for the task at hand. Due to their simplicity, however, 
analytic fits ha ve sometimes been used. As an example, the for- 
mula given in lEisenstein & Hul d 19991) agrees with numerical 
results to better than 5% (around a central region of parame- 
ter spac e). However, using such a fit (as for example, in the 
work of ICrocce et al.ll2010l) gives two errors: (i) an incorrect 
mapping from a(M) to M (Figure IA14I where using the wrong 
transfer function yields a noticeable overprediction of the halo 
abundance at high masses), and (ii) a systematic error in the un- 
derlying simulation which no longer accurately represents the 
correct cosmology. It is therefore important to use the most 
accurate transfer function available to obtain precision results 
from simulations. 

Error due to Finite Force Resolution 

To study the effect of finite force resolution errors quantita- 
tively, we use two simulations with identical initial conditions 
and cosmological parameters. Both the runs have a box size 
of length 1.3 Gpc and 1024 3 particles. One of the simulations 
is run using GADGET-2 with a force resolution of 50 kpc as 
specified in Table [2] The other was run with a particle-mesh 

2 http://camb.info 
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FIG. A14. — Ratio of halo abundances using the fitting formula given 
in Table 3 for two differ ent transfer functions- the analytic fitting form 
of Eisenstein & Hu 1 1999) and that generated using CAMB. The ratio of the 
mass function for two different transfer function choices is shown. At the high 
mass end, the mass function is overestimated by several percent if the analytic 
fit is used. 



(PM) code with a 2048 3 spatial grid, corresponding to a force 
resolution of approximately 700 kpc. We call the two runs G 
and PM respectively. 

Since the runs have identical initial conditions, halos in both 
runs should have approximately the same locations. For each 
halo in run G, we expect t o find a correspondi ng "match" in run 
PM. However, as shown in lLukic et alj d2007l) . for a given force 
resolution, halos below a certain mass (or equivalently contain- 
ing a certain number of particles) are not reliably formed. For 
the PM run, the associated prediction for the minimum number 
of particle s required for a h alo to form is » 400 (see Equa- 
tion 30 in iLukic et al.l 120071) . In order to make sure that the 
mass function is not suppressed due to this effect, the actual 
number of particles in a halo should be larger than this min- 
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imal value. We find that ~ 2000 particles per halo is a good 
choice and hunt for matched halo pairs only above a corre- 
sponding mass of 2 • 10 14 M Q (or halos containing more than 
2500 particles). This being the case, we find that nearly 95% 
of the halos have their centers within ~ 400 kpc (consistent 
with the force resolution of the PM run). The ratio between 
the matched halo masses averaged over all the matched pairs 
is shown as a histogram in Figure lATBI As a function of mass 
bin (lower panel), the peak of the ratio distribution remains es- 
sentially unchanged, with median values ranging from 1.041 to 
1.034. Roughly speaking, the data is therefore consistent with 
an overestimate by approximately 4% in individual halo masses 
in the PM run, independent of the mass. (Note that the mean of 
the ratio distribution has a slightly stronger dependence on halo 
mass, ranging from 1 .06 to 1 .036, resulting from tail effects in 
the distribution.) 

Since the effect is small, a simple linear extrapolation is suf- 
ficient to correct each of the runs in Table [2] in order to predict 
the halo mass in the limiting case. The corrected mass M c is 
given by 

M c /M= 1.0-0.04(e/650kpc). (Al) 

Here M is the uncorrected halo mass and e is the force reso- 
lution measured in kpc of the different runs as specified in Ta- 
ble |2 (We have checked that this formula is consistent with 
another smaller set of PM simulations - 256 3 particles, 1024 3 
mesh, with a force resolution of 334 kpc.) The biggest correc- 
tion needed is « 0.6% for individual halo masses for the A runs 
(with a force resolution of 97 kpc). This results in a systematic 
lowering of approximately 2% in the high-mass tail of the mass 
function (primarily run A) as shown in Figure lA~16l 

Error due to Finite Number of Particles in a Halo 

The reason for the scatter and bias in FOF masses due to the 
finite number of particles in a halo is that the FOF algorithm 
aims to capture the mass of a halo within a certain iso-density 
contour, rendering the halo mass sensitive to the accurate de- 
termination of this boundary. Undersampling of the halo will 
lead to particles on the halo boundary tending to link more to 
particles close by than in the case of a well-sampl ed halo with 
the en d result of overestimating the halo mass (see lLukic et al.l 
|2009j for details). 

In their work on the mass function. lWarren et al. I J2006I) sug- 
gested a correction for the FOF halo mass of the form, rc™ 11 = 
n /i(l _w a°' 6 ), where «/, indicates the number of particles in a 
halo. This adjustment is an empirical finding using a set of 
nested volume simulations. Thi s correction factor has been 
used in recent studies, for example lTinker et al] ( l2008l ); ICrocce et al.l 
d2010l) . In broad ou t line th is result is consistent with the find- 
ings of iLukic et al.l d2009l) . The adjustment formula lowers 
masses for halos with smaller numbers of particles and depends 
only on the number of particles within a halo. Unfortunately, 
this adjustment cannot be applied in all circumstances, as the 
details of the correction can depend on the details of the indi- 
vidual simulations. Additionally, this adjustment should not be 
applied at small particle numbers, where it is known to over- 
compensate. 

The problem with FOF halo masses described above is best 
seen in our simulations by comparing results for the mass func- 
tion in overlapping mass bins across the different-s ized simula- 
tion b oxes (with differing mass resolutions, see also lWarren et al. I 
120061) . This is shown in the upper panel in Figure[2] In our work 
we combine five different box sizes (multiple runs for each box 
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FIG. A 15. — Distribution of the ratio of halo masses in the low resolution 
(PM) and high resolution (G) runs in units of 10 14 Mq. The top panel shows 
the distribution for all halos with M > 1.3 X 10 14 Mq in ran G that have a 
matching halo in the PM ran; more than 95% of the halos have matched pairs 
within 347 kpc. The bottom panel shows the distributions for different mass 
bins demonstrating that the shift in the halo mass between the high and low 
resolution runs is practically independent of mass. 



size) with different mass and force resolution to cover a large 
range of halo masses. The figure shows the ratio of the raw sim- 
ulation results from the five different box sizes with respect to 
our best-fit mass function. In the absence of a systematic bias 
across the boxes, the results should match up within Poisson 
errors, but this is clearly not the case. 

In order to correct for the finite particle sampling of the ha- 
los we investigate different st rategies. If all hal os are NFW, 
then following the analysis in ILukic et all d2009t) we can im- 
plement a correction which takes into account the concentra- 
tion and the number of particles in a halo. The general trend 
is that higher concentration halos should have smaller correc- 
tions while, overall, the corrected halo masses would be lower 
(as in the simulations). Since our force resolution is not suf- 
ficient to reliably determine the concentration of the halos, we 
can use the known mass-concentration relation (including scat- 
ter) to estimate a concentration for each halo. It turns out that 
such a correction lowers the mass function overall by almost 
10% but does not provide an accurate match between differ- 
ent boxes. Therefore, this simple modeling fails to explain the 
quantitative results from the simulations, presumably because 
the concentration estimates are too crude and because of the 
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FIG. A16. — Impact of the co rrec tion for finite force resolution errors on the mass function. The ratio of the simulation results with respect to the best-fit z = 
mass function given in Equation I12t is shown. We display results for different box sizes separately to study possible numerical artifacts in the overlapping regions 
of different boxes. The effect only alters high mass halos (see for a comparison the uncorrected mass function in Figure|2) and leads to a systematic lowering of the 
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FIG . A 17. — Impact of finite sampling errors on the mass function. The ratio of the simulation results with respect to the best-fit z = mass function given in 
Equation U2t is shown. The effect from finite sampling is much larger than the effect due to finite force resolution. 
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irregularities of roughly 20% of the halos. 

We therefore use a conservative strategy of a large minimal 
number of particles per halo, «/, > 400, and apply the empirical 
»/i dependent correc tion. We find that the suggested form in 
IWarren et al. I d2006l) can be slightly improved by changing it 
to: 

„corr „ /i . -0.65\ / a o\ 

n h = n h (l-n h ). (A2) 

The result of applying the correction according to the above 
formula is shown in Figure IAP71 

Finite Volume Effects 

In this appendix we analyze the third finite volume effect 
mentioned in Section[3]in detail, related to the (effective) num- 
ber of independent realizations, i.e., the sample variance. Be- 
cause we are using a simple, independent simulation box tech- 
nique, we compute the variance for each of the runs - separately 
for each fixed value of the simulation volume - over the number 
of realizations of each run. Thus for each run (A, B, C, etc.) we 
calculate 



An 1 



B \ 



£ 



[(n(M,z)-n(M,z)) 2 

N mas 



(A3) 



where N mns is the number of realizations for a particular fixed- 
volume run (e.g., N mns = 10 for the A runs). n is the average 
number density of halos in a bin averaged over all the realiza- 
tions and n is the number density for a single realization. Note 
that this error includes both Poisson fluctuations of halos and 
the separate sample variance contribution. The fractional error 
on /(<x) is then calculated as 

A/a-. = Aita (A4) 

/data "bin 

where we include only Poisson errors and ignore the error in 
determining M\, m due to Poisson fluctuations. 

The fractional errors are shown in Figure IA18I The results 
can be directly compared to si mple estimates o f the v ariance 
using the halo model following iHu & Kravtsovl d2003l) . With 
this approach, one finds 



An 



I 1 b 2 (M,z) 
nV + (2?r) 3 



d*kP(k)W\kR box ), 



(A5) 



where V is the volume of the simulation box, Ri, ox = (3V/(47r)) 1//3 , 
W(x) is the spherical top-hat window function and b(M, z) is the 
large scale bias from Equation Q~8] The first term above is the 
shot noise contribution and the second term gives the sample 
variance contribution to the total variance. Note that the lower 
limit of the integral is specified by the simulation box size. For 
the larger boxes, this is irrelevant, but for the GS runs, the ef- 
fect is significant, as also shown in Figure IaToT which demon- 
strates good agreement between our results and the theoretical 
estimate. Note that by using independent smaller volume runs 
we can reduce the error due to sample variance by a signifi- 
cant amount (of course one should keep in mind that too-small 
boxes will systematically bias the mass function low, a small 
and correctable bias in our case, as shown in Fig. [2]). 
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FIG. A18. — Total fractional error (Poisson plus sample variance) for each 
simulation run. Here, the error bars denote the "error on the errors" which is 
simply the Poisson error. As an example, for the A runs with 10 realizations, 
the "error on the errors" is ~ 30%. The d otted lines represent the correspond- 
ing halo model predictions from Equation lA5l for an infinite volume box, while 
the solid curves represent the same predictions for the actual simulation sizes 
used. The two remain close except for the GS runs; see text for a discussion. 
For all the runs, there is satisfactory agreement between the theoretical esti- 
mates (solid curves) and the actual results from the simulations. 



a E(a) / 3 2 a i E z (a) 



5(l + 8) = 

(Bl) 



Keeping only the linear terms in Equation IB 1 1 gives 



S" + 



3 h &w\ s , 3 n,„ 

a E(a) J 2 a 5 E(a) 2 



5 = 



(B2) 



where E(a) = y/fl m /a 3 + ( 1 - f2 m ) is the Hubble factor at a 
redshift z= 1 /a - 1 . Assuming we want to calculate the value 
of S c at z = 0, we require the sphere to collapse at z = 0. We 
set the initial conditions <5- = 0. We varied S' t and checked that 
the final value of S c is quite insensitive to the value of 5'. We 
solve Equation IB II such that S becomes a large number at z = 
(numerically we attain convergence for S = 10 8 ) for a chosen 
initial perturbation, 8,. We iteratively converge to the value of 
Sj such that collapse happens at z = 0, i.e., 5 = 10 8 at z = 0. 
Once Si is computed by solving Equation IB II we use that value 
to solve Equation IB2I The solution of Equation IB2I gives the 
value of S c at z = 0. Thus, for a given cosmology, i.e., for given 
values of Q, m and w, solving Equations IB 1 1 and |B2 1 gives 8 C at a 
given redshift, z. 



SPHERICAL COLLAPSE CALCULATION OF S^ cdm 

In this section we summarize the calculation of S c for any 
cosmology as predicted by the spherical collapse model. We 
are interested in perturbations in a homogeneous sphere. The 
full nonlinear equation for perturbations ( iPace et al.ll2010t) is 
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